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Abstract 

We study relative dispersion of passive scalar in non-ideal cases, i.e. in situ- 
ations in which asymptotic techniques cannot be applied; typically when the 
characteristic length scale of the Eulerian velocity field is not much smaller 
than the domain size. Of course, in such a situation usual asymptotic quan- 
tities (the diffusion coefficients) do not give any relevant information about 
the transport mechanisms. On the other hand, we shall show that the Finite 
Size Lyapunov Exponent, originally introduced for the predictability problem, 
appears to be rather powerful in approaching the non-asymptotic transport 
properties. This technique is applied in a series of numerical experiments in 
simple flows with chaotic behaviors, in experimental data analysis of drifter 
and to study relative dispersion in fully developed turbulence. 
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It is now well known that the Lagrangian motion of test particles in a fluid 
can be highly non trivial even for simple Eulerian field. Up to now there exist 
powerful methods to study in a rigorous way the asymptotic transport properties 
of passive scalar. On the other hand very often (especially in real world) it is 
not possible to characterize dispersion in terms of asymptotic quantities such 
as average velocity and diff"usion coefficients. This happens, typically, in finite 
domain systems with no large scale-separation between the domain size and 
the largest characteristic Eulerian length; more generally when there is not a 
sharp separation among the characteristic length scales of the system. In this 
perspective, we brieffy review a recently introduced method to approach the 
non-asymptotic properties of transport and mixing. We discuss the relevance 
of the Finite Size Lyapunov Exponent for the characterization of diff"usion. In 
particular we stress its advantages compared with the usual way of looking at 
the relative dispersion at fixed delay time. 



I. INTRODUCTION 

Transport processes play a crucial role in many natural phenomena. Among the many 
examples, we just mention the particle transport in geophysical flows which is of obvious 
interest for atmospheric and oceanic issues. The most natural framework for investigating 
such phenomena is to adopt a Lagrangian viewpoint in which the particles are advected by 
a given Eulerian velocity field u{x,t) according to the differential equation 

^ = uix,t) = v{t), (1) 

where, by definition, v{t) is the Lagrangian particle velocity. 

Despite the apparent simplicity of (|1|), the problem of connecting the Eulerian properties 
of u to the Lagrangian properties of the trajectories x{t) is a very difficult task. In the last 
20 — 30 years the scenario has become even more complex by the recognition of the ubiquity 
of Lagrangian chaos (chaotic advection). Even very simple Eulerian fields can generate very 
complex Lagrangian trajectories which are practically indistinguishable from those obtained 
in a complex, turbulent, fiow ||I|-§1. 

Despite these difiiculties, the study of the relative dispersion of two particles can give 
some insight on the link between Eulerian and Lagrangian properties at different length- 
scales. Indeed, the evolution of the separation R{t) = x^'^\t) — x^^\t) between two tracers 
is given by 

^ = v^^\t) - v^^\t) = u{x^^\t) + R{t),t) -u{x'^^\t),t) (2) 

and thus depends on the velocity difference on scale R. It is obvious from @ that Eulerian 
velocity components of typical scale much larger than R will not contribute to the evolution 
of R. Since, in incompressible ffows, separation R typically grows in time we have the 
nice situation in which from the evolution of the relative separation we can, in principle, 
extract the contributions of all the components of the velocity field. For these reasons, 
in this paper we prefer to study relative dispersion instead of absolute dispersion. For 
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spatially infinite cases, without mean drift there is no difference; for closed basins the relative 
dispersion is, for many aspects, more interesting than the absolute one, which is dominated 
by the sweeping induced by large scale flow. 

There are very few general results on the link between Eulerian and Lagrangian prop- 
erties and only for asymptotic behaviors. Let us suppose that the Eulerian velocity field is 
characterized by two typical length-scales: the (small) scale below which the velocity is 
smooth, and a (large) scale Lq representing the size of the largest structures present in the 
flow. Of course, in most non turbulent flows it will turn out that ~ Lq. 

At very small separations -C we have that the velocity difference in (^ can be 
reasonably approximated by a linear expansion in R, which in most time-dependent flows 
leads to an exponential growth of the separation of initially close particles, a phenomenon 
known as Lagrangian chaos 

(lni?(t)) ^lni?(0) + At (3) 

(the average is taken over many couples with initial separation -R(O)). The coefficient A is the 
Lagrangian Lyapunov exponent of the system [0]. The rigorous definition of the Lyapunov 
exponent imposes to take the two limits -R(O) and then t oo: in physical terms these 
limits amount to the requirement that the separation has not to exceed the scale but for 
very large times. This is a very strict condition, rarely accomplished in real flows, rendering 
often infeasible the experimental observation of the behavior (^. 

On the opposite limit, for very long times and for separations R ^ Lq, the two trajecto- 
ries x^^\t) and x^'^\t) feel two velocities which can practically be considered as uncorrected. 
We thus expect normal diffusion, i.e. 

{R\t))^2Dt. (4) 

Also in this case it is necessary to remark that the asymptotic behavior (^) cannot be attained 
in many realistic situations, the most common of which is the presence of boundaries at a 
scale comparable with Lq. In absence of boundaries it is possible to formulate sufficient 
conditions on the nature of the Eulerian flow, under which normal diffusion (1) always takes 
place asymptotically 0. 

Between the two asymptotic regimes (^ and (^ the behavior of R{t) depends on the 
particular flow. The study of the evolution of the relative dispersion in this crossover regime 
is very interesting and can give an insight on the Eulerian structure of the velocity field. 

To summarize, in all systems in which the characteristic length-scales are not sharply 
separated, it is not possible to describe dispersion in terms of asymptotic quantities. In 
such cases, different approaches are required. Let us mention some examples: the symbolic 
dynamics approach to the sub-diffusive behavior in a stochastic layer [|10[ and to mixing 



in meandering jets [0; the study of tracer dynamics in open flows in terms of chaotic 



scattering [O] and the exit time description for transport in semi-enclosed basins |T3[ and 



open flows [14 



The aim of the present paper is to discuss the use of an indicator - the Finite Size 
Lyapunov Exponents (FSLE), originally introduced in the context of predictability problems 
||T3| - to study and characterize non-asymptotic transport in non-ideal systems, e.g. closed 
basins and systems in which the characteristic length-scales are not sharply separated. 
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In section || we introduce the basic tools for the finite-scale analysis and we discuss their 
general properties. Section III is devoted to the evaluation of our method on some numerical 
examples. We shall see that even in apparently simple situations the use of finite scale 
analysis avoids possible misinterpretation of the results. In section |^ the method is applied 
to two physical problems: the analysis of experimental drifter data and the numerical study 
of relative dispersion in fully developed turbulence. Conclusions are presented in section 0. 
The appendices report, for sake of self-consistency, some technical aspects. 



II. FINITE SIZE DIFFUSION COEFFICIENT 

In order to introduce the finite size analysis for the dispersion problem let us start 
with a simple example. We consider a set of N particle pairs advected by a smooth (e.g. 
spatially periodic) velocity field with characteristic length Z„. Denoting with R'f{t) the square 
separation of the i-th couple, we define 

{R'it)) = ^j:Ri (5) 

We assume that the Lagrangian motion is chaotic, thus we expect the following regimes to 
hold 

\^yt))-\2Dt if (i?2(t))i/2>/^ ' w 



where L{2) > 2A is the generalized Lyapunov exponent |T7,18|, D is the diffusion coefficient 
and we assume that Ri{0) = Rq. 

An alternative method to characterize the dispersion properties is by introducing the 
"doubling time" r(5) at scale 6 as follows ||16|: given a series of thresholds S^""^ = r"5(°) 



one can measure the time Tj(5'^°)) it takes for the separation Ri{t) to grow from S^'^^ to 
(^(1) = r5^'^\ and so on for Ti{5^^^) , Ti(S^'^^) , . . . up to the largest considered scale. The r 
factor may be any value > 1, properly chosen in order to have a good separation between 
the scales of motion, i.e. r should be not too large. Strictly speaking, t{6) is exactly the 
doubling time if r = 2. 

Performing the doubling time experiments over the N particle pairs, one defines the 
average doubling time t{6) at scale 6 as 

1 ^ 

r{s)=<m>^=-j:m). (?) 

It is worth to note that the average (0) is different from the usual time average (see Appendix 



Now we can define the Finite Size Lagrangian Lyapunov Exponent (see [|15| for a detailed 
discussion) in terms of the average doubling time as 

In r 

MS) = ^y (8) 
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which quantifies the average rate of separation between two particles at a distance 5. Let 
us remark that X{S) is independent of r, for r close to 1. For very small separations (i.e. 
6 ^ lu) one recovers the Lagrangian Lyapunov exponent A, 

A = lim — — In r . (9) 

s^ot{6) ^ ' 



In this framework the finite size diffusion coefficient [|I6|, D(5), dimensionally turns out to 
be 

D{6) = 6^X{6) . (10) 

Note the absence of the factor 2, as one may expect from (P), in the denominator of D(S); 
this is because t{6) is a difference of times. For a standard diffusion process D{6) approaches 
the diffusion coefficient D (see eq. (|)) in the limit of very large separations {6 ^ lu)- This 
result stems from the scaling of the doubling times t{5) ~ 5^ for normal diffusion. 
Thus the finite size Lagrangian Lyapunov exponent A (5) behaves as follows: 

^ { )d/5' If 5 >> t ' (11) 

One could naively conclude, matching the behaviors at 5 ~ that D ~ A/^. This is not 
always true, since one can have a rather large range for the crossover due to nontrivial 
correlations which can be present in the Lagrangian dynamics [Q. 

One might wonder that the introduction of t{5) is just another way to look at {R^{t)). 
This is true only in limiting cases, when the different characteristic lengths are well separated 
and intermittency is weak. A similar idea of using times for the computation of the factor 
diffusion coefficient in nontrivial cases was developed in Ref. |n]-pT|. 



If one wants to identify the physical mechanisms acting on a given spatial scale, the use 
of scale dependent quantities is more appropriate than time dependent ones. 

For instance, in presence of strong intermittency (which is indeed a rather usual situation) 
i?^(t) as a function of t can be very different in each realization. Typically one has (see 
figure la), different exponential growth rates for different realizations, producing a rather 
odd behavior of the average {R^{t)) not due to any physical mechanisms. For instance in 
figure lb we show the average {R^{t)) versus time t; at large times one recovers the diffusive 
behavior but at intermediate times appears an "anomalous" diffusive regime which is only 
due to the superposition of exponential and diffusive contributions by different samples at 
the same time. On the other hand, by exploiting the tool of doubling times one has an 
unambiguous result (see figure Ic) [0]. 

An important physical problem where the behavior of r(5) is essentially well understood 
is the relative dispersion in 3D fully developed turbulence. Here the smallest Eulerian scale lu 
is the Kolmogorov scale at which the fiow becomes smooth. In the inertial range lu < R < Lq 
we expect the Richardson law to hold {R^{t)) ~ t^; for separations larger than the integral 
scale Lq we have normal diffusion. In terms of the finite size Lyapunov exponent we thus 
expect three different regimes: 

1. A(5) = A for 5 < 



5 



2. X{S) ~ 5-2/3 for lu<^S<^h 



3. A(5) ~ for 5 > Lo 

We will see in section ^ than even for large Reynolds numbers, the characteristic lengths 
lu and Lq are not sufficiently separated and the different scaling regimes for {R^{t)) cannot 
be well detected. The fixed scale analysis in terms of A (5) for fully developed turbulence 
presents clear advantages with respect to the fixed time approach. 



III. NUMERICAL RESULTS ON SIMPLE FLOWS 

In this section we shall discuss some examples of applications of the above introduced 
indicator A(5) (or equivalently -D(5)) for simple fiows. The technical and numerical details 
of the finite size Lyapunov exponent computation are settled out in Appendix ^ 

In a generic case in addition to the two asymptotic regimes ([TI|) discussed in section 0, 
we expect another universal regime due to the presence of the boundary of given size Lb- 
For separations close to the saturation value 5max — Lb 'v^e expect the following behavior to 



hold for a broad class of systems |16 



^(^) = °^ ^ ■ (12) 

The proportionality constant is given by the second eigenvalue of the Perron-Frobenius 
operator which is related to the typical time of exponential relaxation of tracers' density to 
uniform distribution (see Appendix O). 



A. A model for transport in Rayleigh-Benard convection 

The advection in two dimensional incompressible fiows in absence of molecular diffusion 
is given by Hamiltonian equation of motion where the stream function, ip, plays the role of 
the Hamiltonian: 

dx dip dy dip ,^„s 
dt dy ' dt dx 

If ip is time-dependent one typically has chaotic advection. As an example let us consider 
the time-periodic Rayleigh-Benard convection, which can be described by the following 
stream function |2^: 



iP[x, y,t) = — sin {k[x + B sin(a;t)] }W{y), (14) 

rC 

where W{y) satisfies rigid boundary conditions on the surfaces y = and y = a (we use 
W{y) = sm{7iy/a)). The two surfaces y = a and y = are the top and bottom surfaces of 
the convection cell. The time dependent term B sin(a;t) represents lateral oscillations of the 
roll pattern which mimic the even oscillatory instability [p2| . 

Concerning the analysis in terms of the finite size Lyapunov exponent one has that, if 
5 is much smaller than the domain size, A(5) = A. At larger values of 5 we find standard 
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diffusion X{S) = with good quantitative agreement with the value of the diffusion 

coefficient evaluated by the standard technique, i.e. using as a function of time t. 

In order to study the effects of finite boundaries on the diffusion properties we confine the 
tracers' motion in a closed domain. This can be achieved by slightly modifying the stream 
function (0). We have modulated the oscillating term in such a way that for |a;| = Lb the 
amplitude of the oscillation is zero, i.e. B B sin{7cx/ Lb) with Lb = 27Tn/k [n is the 
number of convective cells). In this way the motion is confined in x G [—Lb, Lb\- 

In figure 2 we show A (5) for two values oi Lb- If is large enough one can distinguish 
the three regimes: exponential, diffusive and the saturation regime eq. (|T^). Decreasing the 
size of the boundary Lb-, the range of the diffusive regime decreases, while for small values 
of Lfi, it disappears. 



B. Point vortices in a Disk 

We now consider a two-dimensional time-dependent ffow generated by M point vortices. 



with circulations Fi, . . . , Fm, in a disk of unit radius The passive tracers are advected 
by the time dependent velocity field generated by the vortices and behave chaotically for 
any M > 2. Let us note that in this case the scale separation is not imposed by hand, but 
depends on M and on the energy of the vortex system p3. Figure 3a shows the relative 



diffusion as a function of time in a system with M = 4 vortices. Apparently there is an 
intermediate regime of anomalous diffusion. However from figure 3b one can clearly see 
that, with the fixed scale analysis, only two regimes survive: exponential and saturation. 
Comparing figure 3a and figure 3b one understands that the appearance of the spurious 
anomalous diffusion regime in the fixed time analysis is due to the mechanism described in 
section 0. 

The absence of the diffusive regime A (5) ~ is due to the fact that the characteristic 
length of the velocity field, which is comparable with the typical distance between two close 
vortices, is not much smaller than the size of the basin. 



C. Random walk on a fractal object: an anomalous diffusive case 

In this section we discuss the case of particles performing a continuous random walk on 
a fractal object of fractal dimension Dp, where one has sub-diffusion. We show that also 
in situation of anomalous diffusion (e.g. sub-diffusion) the FSLE is able to recognize the 
correct behavior. In the following section we consider the case of fully developed turbulence 
which displays super-diffusion. 

In a fractal object due to the presence of voids, i.e. forbidden regions for the particles, 
one expects a decreasing of the spreading, and because of the self-similar structure of the 
domain (i.e. voids on all scales) a sub-diffusive behavior is expected. It is worth to note 
that the particles do not diffuse with the same law from any points of the domain (due to 
the presence of voids), hence in order to define a diffusive-like behavior one has to average 
over all possible particles' position. For discrete random walk on a fractal lattice it is known 
that the diffusion follows the law < R^{t) >~ t^/^w ^j^h Dw > 2, i.e. sub-diffusion 
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25|. The quantity Dw is related to the spectral or fracton dimension Ds, by the relation 



Dw = 2Dp/Ds, and it depends on the detailed structure of the fractal object 



We study the relative dispersion of 2-D continuous random walk in a Sierpinsky Carpet 
with fractal dimension Dp = log 8/ log 3. In our computation we use a resolution 3~^, 
i.e. the fractal is approximated by five steps of the recursive building rule, in practice we 
perform a continuous random walk in a basin obtained with the above approximation of the 
Sierpinsky Carpet. We initialize the particles inside one of the smallest resolved structures, 
then we follow the growth of the relative dispersion with the FSLE method, and redeploy the 
particles in a small cell randomly chosen at the beginning of each doubling time experiment. 
From fig. 4 one can see that X{5) ~ S^'^/-^^ which is an indication of sub-diffusion the 
exponent is in good agreement with the usual relative dispersion analysis (see the inset of 
fig. 4). 



IV. APPLICATION OF THE FSLE 
A. Drifter in the Adriatic Sea: data analysis and modelization 



Lagrangian data recorded within oceanographic programs in the Mediterranean Sea p6 
offer the opportunity to apply the fixed scale analysis to a geophysical problem, for which 
the standard characterization of the dispersion properties gives poor information. 

The Adriatic Sea is a semi-enclosed basin, about 800 by 200 km wide, connected to the 
whole Mediterranean Sea through the Otranto Strait p6| , p7| . We adopt the reference frame 
in which the x,y—axes are aligned, respectively, with the short side (transverse direction), 
orthogonal to the coasts, and the long side (longitudinal direction), along the coasts. 

We have computed relative dispersion along the two axes, (-R^(t)), {Ry(t)) and FSLE 
A (5). The number of selected drifters for the analysis is 37, distributed in 5 different deploy- 
ments in the Strait of Otranto, happened during the period December 1994 - March 1997, 
containing respectively 4, 9, 7, 7 and 10 drifters. To get as high statistics as possible, even 
to the cost of losing information on the seasonal variability, we shift the time tracks of all of 
the 37 drifters to t — Iq, where to is the time of the deployment, so that the drifters can be 
treated as a whole cluster. Moreover, to restrict the analysis only to the Adriatic basin, we 
discarded a drifter as soon as its latitude goes south of 39.5 N or its longitude goes beyond 
19.5 E. 

Before presenting the results of the data analysis, let us introduce a simplified model for 
the Lagrangian tracers motion in the Adriatic Sea. We assume as main features of the surface 



circulation the following elements |2^: the drifter motion is basically two-dimensional; the 
domain is a quasi-closed basin; an anti-clockwise coastal current; two large cyclonic gyres; 
some natural irregularities in the Lagrangian motion induced by small scale structures. On 
the basis of these considerations, we introduce a deterministic chaotic model with mixing 
properties for the Lagrangian drifters. The stream function is given by the sum of three 
terms: 

^(x, y, t) = mo{x, y) + ^i(a;, y, t) + m^ix, y, t) (15) 

defined as follows: 
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c 

^^{x, y) = -^-[- sm{ko{y + vr)) + cos(fco(a; + 2n))] (16) 
C 

"^iix, y,t) = ■ sm{ki{x + sin{ujit))) sm{ki{y + sin(u;it + (/)j))) , [i = 1,2) , (17) 

where ki = 27r/Ai, for 2 = 0, 1,2, Aj's are the wavelengths of the spatial structures of the 
flow; analogously Uj = 27i/Tj, for j = 1,2, and T/s are the periods of the perturbations. 
In the non-dimensional expression of the equations, the length and time units have been set 
to 200 km and 7.5 days, respectively. 

The stationary term \E'o defines the boundary large scale circulation with positive vortic- 
ity. The contribution of \E'i contains the two cyclonic gyres and it is explicitly time- dependent 
through a periodic perturbation. The term gives the motion over scales smaller than the 
size of the large gyres and it is time- dependent as well. The zero-value isoline is defined as 
the boundary of the basin. 

According to observation, we have chosen the parameters so that the velocity range is 
around ~ 0.3 m s~^; the length scales of the Eulerian structures are Lb ^ 1000 km (coastal 
current), Lq ~ 200 km (gyres) and lu ~ 50 km (vortices); the typical recirculation times, 
for gyres and vortices, are ~ 1 month and ~ 1 week, respectively; the oscillation periods are 
~ IQdays (gyres) and ~ 2 days (vortices). 

Let us discuss now the comparison between data and model results. The relative dis- 
persion along the two directions of the basin, for data and model trajectories, are shown in 
figures 5a,b. The results for the model are obtained from the spreading of a cluster of 10^ 
initial conditions. When a particle reaches the boundary (\t' = 0) it is eliminated. 

For the diffusion properties, one cannot expect a scaling for (Rl yit)) before the saturation 
regime, since the Eulerian characteristic lengths are not too small compared with the basin 
size. Indeed, we do not observe a power law behavior neither for the experimental data nor 
for the numerical model. 

Let us stress that by opportunely fitting the parameters, we could obtain the model 
curves even closer to the experimental ones, but this would not be very meaningful since 
there is no clear theoretical expectation in a transient regime. 

Let us now discuss the finite size Lyapunov exponent. The analysis of the experimental 
data has been averaged over the total number of couples out of 37 trajectories, under the 
condition that the evolution of the distance between two drifters is no longer followed when 
any of the two exits the Adriatic basin. 

In fig. 6 we show the FSLE for data and model. In our case, as discussed above, we are 
far from asymptotical conditions, therefore we do not observe the scaling X{S) ~ 



The Xm{S) obtained from the minimal chaotic model (IS-T^) shows the typical step-like 
shape of a system with two characteristic time scales, and offers a scenario about how the 
FSLE of real trajectories may come out. 

The relevant fact is that the large-scale Lagrangian features are well reproduced, at least 
at a qualitative level, by a relatively simple model. We believe that this agreement is not 
due to a particular choice of the model parameters, but rather to the fact that transport is 
mainly dominated by large scales whereas small scale details play a marginal role. 

It is evident the major advantages of FSLE with respect to the usual fixed time statistics 
of relative dispersion: from the relative dispersion analysis of fig. 5 we are unable to recog- 
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nize the underlying Eulerian structures, while the FSLE of fig. 6 suggests the presence of 
structures on different scales and with different characteristic times. In conclusion, the fixed 
scale analysis gives information for discriminating among different models for the Adriatic 
Sea. 



B. Relative dispersion in fully developed turbulence 

We consider now the relative dispersion of particles pairs advected by an incompressible, 
homogeneous, isotropic, fully developed turbulent field. The Eulerian statistics of velocity 
differences is characterized by the Kolmogorov scaling 5v{r) ~ r^/^, in an interval of scales 

-C r -C Lo, called the inertial range, /„ is now the Kolmogorov scale. Due to the 
incompressibility of the velocity field particles will typically diffuse away from each other 
|0J^. For pair separations less than we have exponential growth of the separation of 
trajectories, typical of smooth flows, whereas at separations larger than Lq normal diffusion 
takes place. In the inertial range the average pair separation is not affected neither by large 
scale components of the flow, which simply sweep the pair, nor by small scale ones, whose 
intensity is low and which act incoherently. Accordingly, the separation R{t) feels mainly 
the action of velocity differences 6v{R{t)) at scale R. As a consequence of the Kolmogorov 



scaling the separation grows with the Richardson law [p9| , p0 

(18) 

Non-asymptotic behavior takes place in such systems whenever 1^ is not much smaller 
than Lq, that is when the Reynolds number is not high enough. As a matter of facts even at 
very high Reynolds numbers, the inertial range is still insufficient to observe the scaling (|TH|) 
without any ambiguity. On the other hand, we shall show that FSLE statistics is effective 
already a relatively small Reynolds numbers. 

In order to investigate the problem of relative dispersion at various scale separations 
a practical tool is the use of synthetic turbulent flelds. In fact, by means of stochastic 
processes it is possible to build a velocity fleld which reproduces the statistical properties 
of velocity differences observed in fully developed turbulence |3^. In order to avoid the 



difficulties related to the presence of sweeping in the velocity field, we limit ourselves to 
a correct representation of two-point velocity differences. In this case, if one adopts the 
reference frame in which one of the two tracers is at rest at the origin (the so called a 
Quasi-Lagrangian frame of reference), the motion of the second particle is ruled by the 
velocity difference in this frame of reference, which has the the same single time statistics of 
the Eulerian velocity differences |]32| , |33[| . The detailed construction of the synthetic Quasi- 
Lagrangian velocity field is presented in Appendix |C|. 

In figure 7 we show the results of simulations of pair dispersion by the synthetic turbulent 
field with Kolmogorov scaling of velocity differences at Reynolds number Re ~ 10^ |^|. The 



expected super-diffusive regime (p!8| ) can be well observed only for huge Reynolds numbers 



(see also Ref. |Q. To explain the depletion of scaling range for the relative dispersion let 
us consider a series of pair dispersion experiments, in which a couple of particles is released 
at a separation Rq at time t = 0. At a fixed time t, as customarily is done, we perform an 
average over all different experiments to compute (i?^(t)). But, unless t is large enough that 
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all particle pairs have "forgotten" their initial conditions, the average will be biased. This 
is at the origin of the flattening of for small times, which we can call a crossover 

from initial condition to self similarity. In an analogous fashion there is a crossover for 
large times, of the order of the integral time-scale, since some couples might have reached a 
separation larger than the integral scale, and thus diffuse normally, meanwhile other pairs 
still lie within the inertial range, biasing the average and, again, flattening the curve {R'^{t)). 
This correction to a pure power law is far from being negligible for instance in experimental 
data where the inertial range is generally limited due to the Reynolds number and the 
experimental apparatus. For example, references [^,^| show quite clearly the difficulties 



that may arise in numerical simulations with the standard approach. 

To overcome these difficulties we exploit the approach based on the fixed scale statistics. 
The outstanding advantage of averaging at a fixed separation scale is that it removes all 
crossover effects, since all sampled pairs belong to the inertial range. The expected scaling 
properties of the doubling times is obtained by a simple dimensional argument. The time it 
takes for particle separation to grow from R to 2R can be estimate as T{R) ~ R/6v{R); we 
thus expect for the inverse doubling times the scaling 

/ 1 \ ^ (MM . (19) 

\t{r)/ r ^ ' 

In figure 8 the great enhancement of the scaling range achieved by using the doubling times 
is evident. In addition, by using the FSLE it is possible to study in details the effect of 



Eulerian intermittency on the Lagrangian statistics of relative dispersion. See Ref. |Q for 
a detailed discussion and a comparison with a multifractal scenario. The conclusion that 
can be drawn is that in this case doubling time statistics makes it possible a much better 
estimate of the scaling exponent with respect to the standard - fixed time - statistics. 



V. CONCLUSIONS 



In the study of relative dispersion of Lagrangian tracers one has to tackle situations in 
which the asymptotic behavior is never attained. This may happen in presence of many 
characteristic Eulerian scales or, what is typical of real systems, in presence of boundaries. 



It is worth to stress that such kind of systems are very common in geophysical fiows [|T3l , and 
also in plasma physics [^. Therefore a close understanding of non-asymptotic transport 
properties can give much relevant information about these natural phenomena. 

To face these problems, in recent years, there have been proposed different approaches 
whose common ingredient is basically an "exit time" analysis. We remind the symbolic 
dynamics pi?| , P | and the chaotic scattering |T^ approaches, the exit time description for 
transport in semi-enclosed basins [1^, symplectic maps [0, open flows 0] and in plasma 
physics [PT|]. 

In this paper we have discussed the applications of the Finite Size Lyapunov Exponent, 
A (5), in the analysis of several situations. This method is based on the identiflcation of the 
typical time r(5) characterizing the diffusive process at scale 5 through the exit time. This 
approach is complementary to the traditional one, in which one looks at the average size of 
the clouds of tracers as function of time. For values of 5 much smaller than the smallest 
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characteristic length of the Eulerian velocity field, one has that X{S) coincides with the 
maximum Lagrangian Lyapunov Exponent. For larger 6 the shape of \{S) depends on the 
detailed mechanisms of spreading, i.e. the structure of the advecting velocity field and/or 
the presence of boundaries. The diffusive regime corresponds to the behavior X{5) ~ D/S'^. 
If 5 gets close to its saturation value, i.e. the characteristic size of the basin, the universal 
shape of X{S) can be obtained on the basis of dynamical system theory. In addition, we have 
shown that the fixed scale method is able to recognize the presence of a genuine anomalous 
diffusion. 

A remarkable advantage of working at fixed scale (instead of at fixed time as in the tra- 
ditional approach) is its ability to avoid misleading results, for instance apparent anomalous 
scaling over a certain time interval. Moreover, with the FSLE one obtains the proper scaling 
laws also for a relatively small inertial range for which the standard technique gives rather 
controversial answers. 

The proposed method can be also applied in the analysis of drifter experimental data or 
in numerical model for Lagrangian transport. 
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APPENDIX A: COMPUTATION OF THE FINITE SIZE LYAPUNOV 

EXPONENT 

In this appendix we discuss in detail the method for computing the Finite Size Lya- 
punov Exponent for both continuous dynamics (differential equations) and discrete dynamics 
(maps). 

The practical method for computing the FSLE goes as follows. Defined a given norm for 
the distance S{t) between the reference and perturbed trajectories, one has to define a series 
of thresholds 5„ = r"5o {n = 1, . . . ,P), and to measure the "doubling times" T,.(5„) that a 
perturbation of size 5„ takes to grow up to 5n+i. The threshold rate r should not be taken 
too large, because otherwise the error has to grow through different scales before reaching 
the next threshold. On the other hand, r cannot be too close to one, because otherwise the 
doubling time would be of the order of the time step in the integration. In our examples we 
typically use r = 2 or r = \/2. For simplicity is called "doubling time" even if r 7^ 2. 

The doubling times T,.(5„) are obtained by following the evolution of the separation from 
its initial size 6min Sq up to the largest threshold 6p. This is done by integrating the two 
trajectories of the system starting at an initial distance 5min- In general, one must choose 
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Smin ^ ^0, in order to allow the direction of the initial perturbation to align with the most 
unstable direction in the phase-space. Moreover, one must pay attention to keep 6p < 6max, 
so that all the thresholds can be attained {Smax is the typical distance of two uncorrelated 
trajectory). 

The evolution of the error from the initial value 6min to the largest threshold 6p carries 
out a single error-doubling experiment. At this point one rescales the model trajectory at 
the initial distance 6min with respect to the true trajectory and starts another experiment. 
After error-doubling experiments, we can estimate the expectation value of some quantity 
A as: 

(^)e = ^EA. (Al) 



i=l 



This is not the same as taking the time average because different error doubling experiments 
may takes different times. Indeed we have 

{A). = U" Am = ^^=^-^. (A2) 
T Jo Ei n {r)e 

In the particular case in which A is the doubling time itself we have from 

X{S^) = -J—-\nr. (A3) 

{Ir[On))e 

The method above described assumes that the distance between the two trajectories is 
continuous in time. This is not true for maps or for discrete sampling in time, thus the 
method has to be slightly modified. In this case Tr{6n) is defined as the minimum time at 
which 6(Tr) > r6n- Because now 6(Tr) is a fluctuating quantity, from ( [A^ ) we have 



We conclude by observing that the computation of the FSLE is not more expensive 
than the computation of the Lyapunov exponent by standard algorithm. One has simply to 
integrate two copies of the system and this can be done also for very complex simulations. 



APPENDIX B: UNIVERSAL SATURATION BEHAVIOR OF \{6) 

In this appendix we present the derivation of the asymptotic behavior (|r^) of X{S) for 
S close to the saturation. The computation is explicitly done for the simple case of a one 
dimensional Brownian motion in the domain [—Lb, Lb], with reflecting boundary conditions: 
the numerical simulations indicate that the result is of general applicability. 

The evolution of the probability density p is ruled by the Fokker-Planck equation 

dp_l d^p 

m - 2^d^ ^^^^ 

with the Neumann boundary conditions ^{±Lb) = . 



13 



The general solution of (pl| ) is 

oo 

p{x,t)= ^ 0)e^'=^e-*/"'= + c.c (B2) 

fe=— oo 

where 

rk= (ffr^') ' A; = 0,±l,±2,... (B3) 

At large times p approaches the uniform solution po = ^/'^Lb- Writing p as p{x,t) = 
Po + Sp{x, t) we have, for t ^ ri , 

5p ~ exp(-t/ri). (B4) 
The asymptotic behavior for the relative dispersion {R^{t)) is 

(i?2(t)) = \j[x- x'fp{x, t)p{x\ t) dx dx' (B5) 

For t > Ti using O we obtain {R'^it)) ~ - Ae-^/^^^ Therefore for S^{t) = {R^{t)) 

one has 6{t) ~ — It^^"*^"^^) '^^^ saturation value of 6 is 5max = Lb/V^, so for t ^ ri, 
or equivalently for {6max — S)/6 1, we expect 

(B6) 

which is (|12D. 

Let us remark that in the previous argument for X{S) for S ~ Smax the crucial point is the 
exponential relaxation to the asymptotic uniform distribution. In a generic deterministic 
chaotic system it is not possible to prove this property in a rigorous way. Nevertheless 
one can expect that this request is fulfilled at least in non-pathological cases. In chaotic 
systems the exponential relaxation to asymptotic distribution corresponds to have the second 
eigenvalue a of the Perron-Frobenius operator inside the unitary circle; the relaxation time 
is Ti = — In I a I 



APPENDIX C: SYNTHETIC TURBULENT VELOCITY FIELDS 

The generation of a synthetic turbulent field which reproduces the relevant statistical 
features of fully developed turbulence is not an easy task. Indeed to obtain a physically 
sensible evolution for the velocity field one has to take into account the fact that each eddy 
is subject to the action of all other eddies. Actually the overall effect amounts only to two 
main contributions, namely the sweeping exerted by larger eddies and the shearing due to 
eddies of comparable size. This is indeed a substantial simplification, but nevertheless the 
problem of properly mimicking the effect of sweeping is still unsolved. 

To get rid of these difficulties we shall limit ourselves to the generation of a synthetic 
velocity field in Quasi- Lagrangian (QL) coordinates thus moving to a frame of reference 
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attached to a particle of fluid ri(t). This choice bypasses the problem of sweeping, since 
it allows to work only with relative velocities, unaffected by advection. As a matter of 
fact there is a price to pay for the considerable advantage gained by discarding advection, 
and it is that only the problem of two-particle dispersion can be well managed within this 
framework. The properties of single-particle Lagrangian statistics cannot, on the contrary, 
be consistently treated. 

The QL velocity differences are defined as 

v{r, t)=u (ri(t) +r,t)-u {r^{t),t) , (CI) 

where the reference particle moves according to 

^ = uir,it),t). (C2) 
These velocity differences have the useful property that their single-time statistics are the 



same as the Eulerian ones whenever considering statistically stationary flows . For fully 
developed turbulent flows, in the inertial interval of length scales where both viscosity and 
forcing are negligible, the QL longitudinal velocity differences show the scaling behavior 



r 

v(r) ■ — 

r 



^) ~ r'^" (C3) 



where the exponent is a convex function of and = 1- This scaling behavior is a 
distinctive statistical property of fully developed turbulent flows that we shall reproduce by 
means of a synthetic velocity field. In the QL reference frame the first particle is at rest in 
the origin and the second particle is at r2 = ti + -R, advected with respect to the reference 
particle by the relative velocity 

v{R, t)=u + R,t)-u (ri(t), t) (C4) 

By this change of coordinates the problem of pair dispersion in an Eulerian velocity field 
has been reduced to the problem of single particle dispersion in the velocity difference field 
v{r,t). This yields a substantial simplification: it is indeed sufficient to build a velocity 
difference field with proper scaling features in the radial direction only, that is along the 
line that joins the reference particle ri{t) - at rest in the origin of the QL coordinates - 
to the second particle r2(t) = ri{t) + R{t). To appreciate this simplification, it must be 
noted that actually all moments of velocity differences u {ri{t) + r',t) — u {ri{t) + r,t) = 
v{r', t) — v{r, t) should display power law scaling in |r' — r | . Actually these latter differences 
never appear in the dynamics of pair separation, and so we can limit ourselves to fulfill the 
weaker request Needless to say, already for three particle dispersion one needs a field 

with proper scaling in all directions. 

We limit ourselves to the two-dimensional case, where we can introduce a stream function 
for the QL velocity differences 

i;(r,t) = V X ^(r,t) . (C5) 

The extension to a three dimensional velocity field is not difficult but more expensive in 
terms of numerical resources. 
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Under isotropic conditions, the stream function can be decomposed in radial octaves as 



^(r, ^, t) = E E ^^F{hr)G,,{e) . (C6) 

i=l j=l '^i 

where ki = 2\ Following a heuristic argument, one expects that at a given r the stream 
function is essentially dominated by the contribution from the i term such that r ~ 2~\ 
This locality of contributions suggests a simple choice for the functional dependencies of the 
"basis functions" : 

F{x) = x^{l -x) for < a; < 1 (C7) 



and zero otherwise. 



G^,l{0) = 1, G,,2{0) = cos(2e + <^,) (C8) 



and Gij = for j > 2 {ipi is a quenched random phase). It is worth remarking that this 
choice is rather general because it can be derived from the lowest order expansion for small 
r of a generic streamfunction in Quasi-Lagrangian coordinates. 

It is easy to show that, under the usual locality conditions for infra red convergence, (p < 
p [^], the leading contribution to the p-th order longitudinal structure function {\vr{r)\P) 
stems from M-th term in the sum (|C6|), {\vr{r)\^) ~ {\4>m,2\^) with r ~ 2~*^. If the (f)ij{t) are 
stochastic processes with characteristic times x; = 2"^*/'^ro, zero mean and {\4>i,j\^) ~ fcj 



the scaling (JC^) will be accomplished. An efficient way of to generate is |^ 



Mi) = 9^At) ^2,,(t) ■ ■ ■ Z,,j{t) (C9) 

where the Zkj are independent, positive definite, identically distributed random processes 
with characteristic time r^, while the gi^j are independent stochastic processes with zero 
mean, {gf^) ~ and characteristic time Tj. The scaling exponents Cp are determined by 
the probability distribution of via 

Cp = ^-log2(z0. (CIO) 

As a last remark we note that by simply fixing the = 1 we recover the Kolmogorov 
scaling, which has been used in the simulations presented in section |1V B 
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FIGURE CAPTIONS 



Figure 1: a) Three realizations of -R^(t) as a function of t built as follows: R'^{t) = 
5^exp(27t) if R^{t) < 1 and R'^{t) = 2D{t - Q with 7 = 0.08,0.05,0.3 and 
60 — 10~^, D — 1.5. b) (-R^(i)) as function of t averaged on the three realizations 
shown in figure la. The apparent anomalous regime and the diffusive one are shown, 
c) X{S) vs S, with Lyapunov and diffusive regimes. 

Figure 2: Lagrangian motion given by the Rayleigh-Benard convection model with: A — 
0.2, B = 0.4, uj = 0.4, A; = 1.0, a = tt, the number of realizations is A/" = 2000 
and the series of thresholds is 5„ = 5o^" with So = 10^'^ and r = 1.05. \{6) vs 6, in 
a closed domain with 6 (crosses) and 12 (diamonds) convective cells. The lines are 
respectively: (a) Lyapunov regime with A = 0.017; (b) diffusive regime with D — 0.021; 
(c) saturation regime with Smax — 19.7; (d) saturation regime with Smax — 5.7. 

Figure 3: (a)(i?^(t)) for the four vortex system with Fi = r2 = — Fs = — F4 = 1. The 
threshold parameter is r = 1.03 and Sq — 10~^, the dashed line is the power law 
{R'^{t)) ~ t^-^. The number of realizations is Af = 2000. (b)A(5) vs 5 for the same 
model and parameters. The horizontal line indicates the Lyapunov exponent (A = 
0.14), the dashed curve is the saturation regime with Smax = 0.76 

Figure 4: FSLE computed for particle diffusion in a Sierpinsky Carpet of fractal dimension 
Df = log(8)/ log(3) obtained by iteration of the unit structure up to a resolution 3~^, 
one has: X{6) ~ 6^^^'^^, which is in agreement with the value obtained for < R{t) > 
versus t shown in the inset (i.e. < R{t) >~ t'^^). 

Figure 5: Relative dispersion of Lagrangian trajectories, < R'^y{t) > versus t, in the Adri- 
atic Sea, for data (continuous line) and model (dashed line), along the natural axes of 
the basin: (a) transverse direction (x— axis) and (b) longitudinal direction (y— axis). 
The time is measured in days and the mean square radius of the cluster is in krn^. 

Figure 6: FSLE of Lagrangian trajectories in the Adriatic Sea, for data (continuous line) 
and model (dashed line). The scale 5 is in km, \{S) is in days"^. 

Figure 7: Relative dispersion {R{t)) for iV = 20 octaves synthetic turbulent simulation 
averaged over 10^ realizations. The line is the theoretical Richardson scaling t^/^. 

Figure 8: Average inverse doubling time (1/T(i?)) for the same simulation of the previous 
figure. Observe the enhanced scaling region. The line is the theoretical Richardson 
scaling 
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